Activated volcanism of Mount Fuji by the 2011 Japanese large earthquakes

The relation between earthquakes and volcanic eruptions, each of which is manifested by large-scale tectonic plate and mantle motions, has been widely discussed. Mount Fuji, in Japan, last erupted in 1707, paired with a magnitude (M)-9-class earthquake 49 days prior. Motivated by this pairing, previous studies investigated its effect on Mount Fuji after both the 2011 M9 Tohoku megaquake and a triggered M5.9 Shizuoka earthquake 4 days later at the foot of the volcano, but reported no potential to erupt. More than 300 years have already passed since the 1707 eruption, and even though consequences to society caused by the next eruption are already being considered, the implications for future volcanism remain uncertain. This study shows how volcanic low-frequency earthquakes (LFEs) in the deep part of the volcano revealed unrecognized activation after the Shizuoka earthquake. Our analyses also show that despite an increase in the rate of occurrence of LFEs, these did not return to pre-earthquake levels, indicating a change in the magma system. Our results demonstrate that the volcanism of Mount Fuji was reactivated by the Shizuoka earthquake, implying that this volcano is sufficiently sensitive to external events that are considered to be enough to trigger eruptions.


Results
In this section, we examine the temporal change in the activity of LFEs beneath Mount Fuji, Japan, and discuss the possibility of an increase in magmatic activities in the lower crust which was triggered by the 2011 Tohoku and Shizuoka earthquakes. We first obtain new catalogs of LFEs, and then apply statistical analyses to examine the temporal changes in LFE activity. Based on the results of these analyses, we conclude that the change in stress caused by the 2011 Shizuoka earthquake enhanced the creation of fractures and resultant change in activity of LFEs around the deep magma reservoir in the lower crust. We also suggest that the activation of LFEs may reflect the sensitivity of Mount Fuji to external disturbances and that the volcano may have the potential to erupt after a 300-year quiescent period.
First, we show the influence of the 2011 Japanese earthquakes on LFEs manifested in the magma system of Mount Fuji. Second, we offer support that this magma system was disturbed due to such events, which played a role on activated LFEs. For inference on volcanic hazards of Mount Fuji, implications for future volcanism based on our observations will be discussed later.
Influence on LFE activity after the 2011 Japanese earthquakes. To resolve the difficulty in detecting LFEs by conventional event-detection methods, we produced seismicity catalogs using the matched-filter (MF) method, which cross-correlated a template to continuous seismic signals [22][23][24][25][26] (details of the "MF method" in "Methods" and Supplementary Figs. S1-S4). In this method, which considers a continuum of seismic signals, an LFE is identified when the timing of a seismic signal and the timing of a template signal overlap. Our resultant catalog includes ~ 16,000 LFEs ( Supplementary Fig. S1) whose correlation coefficients (CC) to quantify signal similarity range from 0.1 (poor correlation) to 1 (strong correlation and most likely a template LFE). Histograms of CC-values show a tall peak at CC ~ 0.2 ( Supplementary Fig. S1). The minimum threshold for CC (CC th ), above which LFEs are used for subsequent analysis, should be above the upper noise limit. A low CC th would result in a large number of false detections while a high CC th would result in a reduced number of available LFEs. We selected CC th = 0.25, at the 2-sigma level of a normal distribution, as is expected for the random correlation between signal and noise [27][28][29] . Higher values of CC th = 0.3, 0.35, and 0.5 above the 3-sigma level were also considered to examine the dependence of the result on CC th . Mainly the results for CC th = 0.25 are shown in the text.
The number of LFEs in our catalog of CC th = 0.25 is about three times larger than in the JMA catalog, which lists about 2,000 LFEs detected in 2003-2019 by a conventional method that is not based on CC (Figs. 1c and 2a). Our catalogs of CC th = 0.25, 0.3, and 0.35 include LFEs immediately after the Shizuoka earthquake ( Fig. 2a and Supplementary Fig. S1) while the JMA catalog does not (Fig. 1d). One may see that there are just a handful of LFEs next to the line that indicates the timing of the Shizuoka earthquake (Fig. 2b). There may be doubts as to whether high-frequency (ordinary) aftershocks were mis-detected as LFEs, since the source area of the Shizuoka earthquake and the occurrence area of LFEs were close to each other (Fig. 1b). If CC th was selected to be below the upper noise limit, there would be a chance that our catalogs inadvertently included many non-volcanic seismic events. However, our choice of CC th (= 0.25, 0.3, 0.35) indicates that events with low-frequency signals truly occurred immediately after the Shizuoka earthquake. Given that the number of LFEs in our CC th = 0.5 catalog is about one-fourth of that in the JMA catalog ( Supplementary Fig. S1), there were no LFEs immediately after the Shizuoka earthquake.
The cumulative number of LFEs as a function of time in Fig. 2a shows that the rates of occurrence changed around the moments of the Tohoku and Shizuoka earthquakes. To support the observation that either earthquake very likely played a role in reactivation, a statistical analysis was conducted. We introduced the epidemic type aftershock sequence (ETAS) model 30 , assuming that this model, which was originally developed for ordinary earthquakes, was applicable to LFEs (details of the ETAS model in "Methods"). Even though the ETAS model provides a good fit to standard earthquake occurrence, it is known that transient non-standard cases are poorly fitted by the standard ETAS model [31][32][33] . We focused on differences between the standard ETAS model and an extended two-stage ETAS model, which covers non-standard cases in fitting LFE timeseries, whereas the twostage ETAS model is the simplest alternative to the standard model 32,33 . This two-stage model assumes different parameter values in subperiods before and after a particular time (change point, T c ), and the whole period is www.nature.com/scientificreports/ divided into two adjoining periods to fit the ETAS models separately. This is best applied to cases where there is a clear-cut time instant across which ETAS parameter values change 32,33 . To test whether or not there are changes in seismicity pattern at T c in a given period, the problem of model selection is reduced by using AIC (Akaike Information Criterion) 34 . We compared AIC single (AIC for a single ETAS fitting) with AIC 2stage (AIC for a twostage ETAS fitting) to select the model with the smaller value, where AIC 2stage = AIC 1 + AIC 2 (AIC 1 and AIC 2 for fitting to the 1st and 2nd subperiods, respectively). ΔAIC (= AIC single − AIC 2stage ), as a function of T c during 2003-2019 (Fig. 2b), shows that the two-stage ETAS model was much better than the single ETAS model, indicating that the most significant T c was around the time of the Tohoku and Shizuoka earthquakes. This feature was observed regardless of the CC th value ( Supplementary  Fig. S5). There was a pronounced discontinuation in smoothing at the time of the Shizuoka earthquake due to triggered LFEs for CC th = 0.25 (Fig. 2b), 0.3, and 0.35, but this was not the case for CC th = 0.5 due to few LFEs around the time of the Shizuoka earthquake. AIC when the JMA catalog was used (grey data in Fig. 2b) also showed a better (but insignificant) outcome when the two-stage ETAS model, rather than the single ETAS model, was used ( Supplementary Fig. S5). Namely, ΔAIC was higher than 0 (ΔAIC > 0) when T c was around the time of the Shizuoka earthquake, but it was below the horizontal dashed line in Fig. 2b. This line indicates a hurdle to the selection of the two-stage ETAS model, given that T c was searched from the data during 2003-2019 32,33 .
Visual inspection of T c immediately before the Shizuoka earthquake confirmed that fitting by the two-stage ETAS model (Fig. 3a,b) was better than that by the single ETAS model (Fig. 3c). An AIC 1 of 216.7, when seismicity since 2003 until immediately before the Shizuoka earthquake was fitted, shows that the occurrence rates (black) were significantly larger than the extrapolated rates (red) after the Shizuoka earthquake (Fig. 3a). This significant difference continued until more recent times (until 2019). This continuation (AIC 2 = − 4234.6) (Fig. 3b)  is that the magma system, which is underground, i.e., beneath Mount Fuji, and directly unobservable, results in highly clustered and heterogeneous occurrence times of LFEs. A solution to this problem is to take a reductionism approach, by decomposing LFEs into primary LFEs and secondary-triggered LFEs. Measuring the occurrence rate of the former LFEs (indicative of Poisson background activity) and the latter LFEs (indicative of aftershock activity), allows the magma system and the interactions (triggering) among LFEs to be inferred. The ETAS model can be applied to offer a solution [30][31][32][33] . In the ETAS model 30 , seismic activity is expressed by two terms, one for Poisson background activity and another for aftershock activity, assuming that each earthquake (including the aftershocks of another earthquake) is followed by aftershocks (details of the ETAS model in "Methods").
As a preliminary analysis, we examined whether decomposing the LFE sequence was really meaningful, i.e., if secondary triggering played a role in the LFE sequence. Using AIC for the two-stage ETAS model (same model as for Fig. 2b and Supplementary Fig. S5) and the two-stage Poisson model (same as the two-stage ETAS model except that the aftershock activity term was ignored), we compared the goodness-of-fits between them applied to the same dataset. For this comparison, the same procedure as for Fig. 2b and Supplementary Fig. S5 was used except that the single ETAS model was replaced by the two-stage Poisson model. The difference in AIC between the ETAS and Poisson models shows that the former model is superior to the latter one ( Supplementary Fig. S6), indicating that the aftershock activity term is not negligible. This result allowed us to examine whether or not the changes in occurrence rate were due to enhanced Poisson background and aftershock activities ( Fig. 4 and Supplementary Table S1).
We used μ (a solo parameter of the Poisson background activity term of the ETAS model) and K 0 (a parameter representing clustering aftershock productivity and one of the four parameters of the aftershock activity term), while the other three parameters were constants. For details of these parameters, see the ETAS model in "Methods". In Fig. 4, μ and K 0 were significantly larger after than before 2011 (details of time-dependent μ and K 0 in "Methods" and Supplementary Table S1). Moreover, μ and K 0 increased with decreasing CC th for both periods, before and after 2011. This indicates that small-CC LFEs identified by relaxing the event identification protocol contributed to both the Poisson background activity and the aftershock activity. The μ-pattern was not sensitive to switching from K 0 ( Fig. 4 and Supplementary Table S1) to other parameters (details of time-dependent μ and K 0 in "Methods" and Supplementary Fig. S7 and Tables S2 and S3). Our results show the contribution of enhanced background and aftershock activities to the increase in occurrence rate of LFEs, and support the possibility of changes in the magma system due to the Shizuoka earthquake.
The bottom panel of Fig. 2a shows that the high rate of LFEs after the Shizuoka event decreased as a function of time, implying that LFEs behaved as aftershocks of the Shizuoka event. Therefore, it may be considered that the background rate of LFEs was globally similar to what was observed before 2011. However, this is not true because μ before 2011 was smaller than after 2011 (Fig. 4a) and because the two-stage ETAS model performed better than the single ETAS model, which included the Shizuoka earthquake (Fig. 2b). A clear change in LFE activity during the occurrence of the Shizuoka earthquake indicates the contribution of the changes in both aftershock and background activities.
If there was to be a disturbance in the magma system, it might induce differences in template LFEs between pre-and post-Shizuoka-quake sequences, allowing differences in detected LFEs between them. A simple test (top panel of Fig. 5) was conducted to assess this possibility, revealing that templates in the pre-(post-) Shizuokaquake sequence were likely applicable for detecting LFEs in the same sequence. Figure 5 Supplementary Fig. S8). LFEs before/after 2011 were mostly detected by templates in the same time periods. In this test (top panel of Fig. 5), the two time-windows before and after the Shizuoka earthquake were chosen. It is quite usual for detections to be temporally clustered with template LFEs, so there may be curiosity about the time-windows that were selected. However, the statement that LFEs before/after 2011 were mostly detected by templates in the same time periods is insensitive to selection bias of the time-windows ( Supplementary Fig. S9).    www.nature.com/scientificreports/ volatile components (H 2 O and CO 2 ) from liquid to gas [37][38][39][40][41] . If the magma plumbing system is open or has weak walls that are easily expanded or fractured, increased pressure creates paths for the magma to migrate, so magma pressure will be further reduced and bubbling will be accelerated. In the second process, cracks close to the magma system would be unstable due to stress perturbation. These failures can become paths through which magma flows up, leading to an eruption. Disturbances caused by external fractures, such as the Shizuoka earthquake, may be a key to initiating both processes. Nakamichi et al. 35 suggested that a complex process occurred beneath Mount Fuji in which the characteristics of its LFEs were due to a variety of focal processes, and where the source mechanism of the largest LFE (M2.3) was explained by non-double-couple components. Among variable LFEs, only those associated with newly-created fractures and reactivated fractures pushed closer to failure by stress changes due to the Shizuoka earthquake, may have become active. The reader might note alternative source processes. LFEs might not always be generated by magma movement. Hydrothermal fluids, as well as the attenuation of higher frequency earthquakes, can create apparent LFEs 42,43 .
Schematic cross-sections beneath Mount Fuji 17, 44 before and after the Shizuoka earthquake (Fig. 6) show how newly created fractures 19 and activated LFEs are associated with changes in the magma system and fluid-and gas-rich area 35,36 . However, based on our result, it remains unsolved if the Shizuoka earthquake triggered an increase in the supply of magmatic fluid to the magma reservoir. Furthermore, there was no evidence that the amount of rising magma increased because clear near-surface deformation around Mount Fuji was not observed before and after the Shizuoka earthquake.
Stress changes on the magma system beneath Mount Fuji of 0.001-0.01 and 0.1-1 MPa were caused by the Tohoku and Shizuoka earthquakes, respectively 19 . The latter change is considered to be sufficient to trigger earthquakes 20,21 , implying the excitement of the magma system and triggering an eruption through either or both processes described above 19 . Our results (Figs. 2, 3, 4, 5) demonstrate that the Shizuoka earthquake played a role in the magma system's excitation, but not enough to trigger an eruption. We conclude that Mount Fuji was sensitive to disturbances due to this earthquake. This is consistent with a previous study in which it was implied that the crust in the area near Mount Fuji is quite sensitive to transient stress perturbation and that the level of pressurization of the hydrothermal and/or magmatic fluids is high in the Mount Fuji area 45 .
In 1703, four years before the 1707 Hoei eruption, seismic swarms were observed 35 days after the M8-class Genroku Kanto earthquake, about 100 km east of Mount Fuji, but no eruption occurred [14][15][16] . The 1707 eruption was also preceded by the M9-class Hoei earthquake, about 200 km to the southwest, on October 7, 49 days before the eruption. Beginning on November 28, 1707, earthquake swarms were observed several times and dozens of earthquakes were felt from December 15, 1707. Mount Fuji then began to erupt on December 16, 1707.
In 2000-2001, LFE swarms occurred, starting in August 2000, two months after the eruption of Miyakejima volcano, which lies 160 km to the south of Mount Fuji, although the change in stress was 10 -4 MPa 19, 46 , or 0.001 ~ 0.0001 of that caused by the Shizuoka earthquake. This change in stress is considered to be too small to trigger an eruption. The experience of Mount Fuji described above implies that it tends to be influenced by external disturbances such as large earthquakes and active volcanoes. Our observation of activated LFEs due to the Shizuoka earthquake and a previous theory of an increase in stress imparted by this earthquake 19 support this tendency, although the volcano has not yet erupted (June 2023). While this study presents a case for the interaction of the Mount Fuji volcanic system with tectonic earthquakes, it remains possible that considerable volcano-seismic activity took place without any influence by tectonics for some cases.
Over 300 years after the 1707 Hoei eruption, the Japanese government has started to consider preparations for the next eruption, citing a worst-case scenario that inflicts catastrophic damage on humans and society 18 . Whether the increase in occurrence rate of LFEs (Fig. 2) will continue or subside absent external events remains unknown, and it is impossible to draw conclusions about the timing of the next eruption. However, our detailed study demonstrates that LFE activity is an important indicator of Mount Fuji's subsurface magma system 17,35,36,44 . Given that this study covered data up to 2019, additional analyses for more recent LFEs in future research

MF method.
When studying LFEs associated with volcanic phenomena, researchers may want to use a catalog that consists of a complete list of LFEs. However, due to their low signal-to-noise ratio, LFE signals are difficult to detect by conventional event-detection methods. We used an MF method that correlates waveforms of continuous signals with those of a template and allows the detection of seismic sequences with a low signal-tonoise ratio [22][23][24][25][26] . In this study, the MF system used for detecting LFEs beneath the Hakone volcano, Japan 22 was modified so that it was applicable to Mount Fuji. Waveforms of continuous signals that were used in this study covered the Jan. 2003-Jul. 2019 period, as recorded by 16 seismic stations (Supplementary Fig. S1) with a three-component velocity seismometer around Mount Fuji. These were obtained from the Earthquake Research Institute at the University of Tokyo.
All stations used in this study ( Supplementary Fig. S1) host borehole seismometers (eigen frequency of 1 Hz), except for the stations OMZ To prepare template LFEs, we used the JMA catalog, which includes ordinary earthquakes and LFEs observed in Japan. Although ordinary earthquakes are distributed all over Japan, LFEs tend to concentrate beneath active volcanoes (Fig. 1a), along the boundary between the Philippine Sea plate and the continental plate in western Japan 47 , and as several isolated clusters in the intraplate regions 48 . Each event in the JMA catalog is classified based on subsidiary information: natural (ordinary) earthquake, LFE, artificial event, etc. The spatial map of events classified as LFEs shows that the area of LFEs around Mount Fuji was separated from other areas of LFEs (Fig. 1a). We only selected LFEs around Mount Fuji and defined the catalog including these LFEs (Fig. 1b-d).
The source area of the Shizuoka earthquake and the occurrence area of LFEs are close to each other (Fig. 1b). However, due to subsidiary information, there is no doubt that the LFE catalog we used eliminated aftershocks (ordinary earthquakes) of the Shizuoka earthquake.
This study relied on statistical analyses of the LFE catalog, which required the use of a complete LFE catalog that covered the study region and time period. It should be carefully considered that the catalog may be controlled by the selection of template earthquakes in the MF analysis. To ensure the completeness of the catalog, it is critical to use a well-chosen set of template LFEs. Careful consideration was needed to make a set of template LFEs, as explained below.
First, large LFEs were selected as templates in order to allow template waveforms to include more information on signals than on noise. Using the JMA catalog, we investigated an M-time graph of LFEs around Mount Fuji during the Jan. 2003-Jul. 2019 period (Fig. 1c), and found that the majority was M = 0 ~ 1, regardless of the date. Visual inspection shows that LFEs of M > 1 were rare with no particular tendency such as LFEs of M > 1 which occurred more frequently or less frequently over time. If smaller magnitude criteria were selected to increase the number of template LFEs, then more LFEs would be detected. However, since implementation of the MF system was computationally more intensive when using a large number of template LFEs than when using a small number of template LFEs, so an implementation trial was conducted by using different magnitude criteria under our computing environment. We found that the most feasible was to use M ≥ 0.9 LFEs as templates. Thus, irrespective of time, LFEs with M ≥ 0.9 in Jan. 2003-Jul. 2019 were selected as templates. Then, among them, LFEs that were recorded by six stations with a minimum signal-to-noise-ratio of 2, were selected 22 .
Second, we verified whether a chosen-set of LFEs (M ≥ 0.9) showed appropriate spatiotemporal coverage. This becomes particularly important when the source mechanism changes over time. Map and cross-sectional views around Mount Fuji (Fig. 1b) show that selected LFEs (M ≥ 0.9) were mostly distributed within the cluster of LFEs, with all magnitudes between Jan. 2003 and Jul. 2019. A consistent spatial coverage of LFEs (M ≥ 0.9) was found between the pre-and post-Shizuoka-quake periods: the coverage of LFEs (M ≥ 0.9) before the Shizuoka earthquake was observed in the area where LFEs (M ≥ 0.9) had already occurred thus far (Supplementary Fig. S1).
It appears that there were more template LFEs from the post-Shizuoka-quake period than from the pre-Shizuoka-quake period ( Fig. 5 and Supplementary Fig. S8), which could potentially result in an overrepresentation of detected earthquakes in the post-Shizuoka-quake period relative to the pre-Shizuoka-quake period. However, as described above, this bias was not intentionally included. Thus, we did not distort to select template LFEs for making a template LFE catalog, which would be used for a subsequent MF procedure.
The MF procedure to identify LFEs, briefly described in this paragraph, is the same as that of Yukutake et al. 22 . Three-component waveform records for each template LFE were used, applying a six-second time window beginning two seconds before the onset time of the theoretical S-wave arrivals. Both templates and continuous waveforms were bandpass-filtered for 1-6 Hz and decimated at 20 Hz to reduce the calculation cost. This band was selected because Yukutake et al. 22 was followed, although other studies used a slightly narrower band (e.g., [1][2][3][4]  www.nature.com/scientificreports/ enhance the reliability of this study. After removing multiple counts, the location of the candidate was assigned to the hypocenter of the matched template LFE determined by JMA (also see the "Catalog quality evaluation" section in "Methods"). Magnitude was determined as the mean of the maximum amplitude ratios of the template with respect to the candidate. The MF procedure described above was applied to all waveform records in the Jan. 2003-Jul. 2019 period, and a preliminary catalog, including candidate LFEs, was created, but LFEs identified by five or less stations were not included in this catalog 22 . Less reliable LFEs were removed from this preliminary catalog to create a finalized catalog, as follows. Among candidate LFEs, false detection occasionally occurred due to contamination by other seismic signals such as teleseismic earthquakes. This contamination led to the detection of LFEs with a large M, so we visually inspected whether each template LFE was used to detect many candidate LFEs with M > 2, a magnitude above which few LFEs have been recorded beneath Mount Fuji in the JMA catalog since 2003. We considered that such template LFEs had a feature similar to teleseismic earthquakes and decided to eliminate them from the list of template LFEs. Thus, candidate LFEs detected by using the eliminated template LFEs were removed from the preliminary catalog, resulting in the finalized catalog that included ~ 16,000 LFEs. A total of 87 template LFEs were used for the finalized catalog. The locations of template LFEs and seismic stations are indicated in Fig. 1b and Supplementary Fig. S1. Despite this primary quality test, an additional test was conducted, as described in the next paragraphs and in the "Catalog quality evaluation" section.
The CC-values of LFEs ( Supplementary Fig. S1) ranged between 0.1 (poor correlation with a template LFE) and 1 (strong correlation with, and identical to, the corresponding template LFE). Setting the minimum CC to a low value implies the use of an incomplete catalog influenced by the nature of low signal-to-noise ratios of LFEs. The minimum threshold for CC (CC th ), above which LFEs are used for our analysis, should be above the upper noise limit. We decided to use CC th = 0.25, 0.3, 0.35, and 0.5 for the following reasons. Histograms of CC-values in Supplementary Fig. S1 show an asymmetric distribution with a tall peak at CC ~ 0.2. We followed previous studies [27][28][29] , in which the distribution of lower CC-values was modeled by a normally distributed curve that would be expected for random correlations between signals and noise, while the upper tail was considered to represent the presence of well-correlated LFEs. Visual inspection shows that frequencies at and below CC ~ 0.2 are in good agreement with the left-hand side of the normally distributed curve where the mean is 0.19 and its standard deviation is 0.03 ( Supplementary Fig. S1). We selected CC th = 0.25, which corresponds to the mean plus two standard deviations. Moreover, we selected CC th = 0.3 and 0.35, which are higher than the mean plus three standard deviations, to examine the dependence of the result on CC th . Similar to Green and Neuberg 27 and Petersen 28 , we found an outliner peak at CC ~ 0.3 ( Supplementary Fig. S1). This peak was clearly observed in the histogram of CC-values for M ≥ 0 ( Supplementary Fig. S1). This histogram is displayed because our analysis basically did not include LFEs with M < 0. Previous researchers, who studied Shishaldin volcanoes (Alaska), the Soufriere volcano (West Indies), and the Unzen volcano (Japan) [27][28][29] , selected CC th -values > 0.5, by showing normally distributed curves with larger means and standard deviations than those shown in this study. We also examined the case for CC th = 0.5.
The scope of this study did not permit us to reveal repeating LFEs, nor cyclic activities and cluster characteristic, as were studied by Lamb et al. 29 . Rather, this study's objective was to resolve the difficulty in detecting smaller LFEs. Our future research will be to conduct in-depth analyses of repeating LFEs for each cluster beneath Mount Fuji, referring to Lamb et al. 29 , and using a sophisticated MF method that can locate detected LFEs to appreciate if they occurred in the same cluster as the template LFE used to find them. Catalog quality evaluation. As a basis of catalog completeness, understanding magnitude scales used in this study is critical. We examined whether large LFEs in our catalog were indeed large, as in the JMA catalog, and vice versa. An LFE in our catalog was paired with that in the JMA catalog if the time difference between them was within two seconds, while ignoring one-to-multiple cases. In this pairing, differences in the locations of LFEs were not considered because, as described in the "MF method" section in "Methods", the locations of LFEs in our catalog were assigned to the hypocenter of the matched template LFE determined by JMA. A list of paired LFEs shows that magnitude in our catalog was positively correlated with that in the JMA catalog, and that the former was nearly equal to the latter ( Supplementary Fig. S4), allowing us to assume a one-to-one transformation in magnitude between our catalog and the JMA catalog.
Analyses of the ETAS model (see the "ETAS model" section in "Methods") are critically dependent on a robust estimate of completeness magnitude (M c ) of the processed LFE data. Above M c , all events are considered to be detected. In particular, underestimates of M c lead to unreliable ETAS fitting. Attention always needs to be paid to M c when assessing M c in each time window. Details of how to compute M c are provided in the "Computation of M c " section in "Methods".
M c was about 0.3 ~ 0.5 for the JMA catalog and about 0.2 ~ 0.3 for the MF catalogs ( Supplementary Fig. S4). These estimates of M c were based on precut catalogs covering several time periods (Supplementary Fig. S4). Therefore, we did not consider a single M c over the entire catalog. To verify whether the results depended on the choice of minimum magnitude (M th ), above which the ETAS model was fitted, we assumed M th = 0. www.nature.com/scientificreports/ izes seismic activity or earthquake productivity of a region, and constant b is used to describe the relative occurrence of large and small events (i.e., a high b-value indicates a larger proportion of small earthquakes, and vice versa). Changes in b-values of ordinary earthquakes are known to reflect structural heterogeneity, strength, and temperature in the Earth [51][52][53][54][55] , and the b-value is also known to be inversely dependent on differential stress 56,57 . We assumed that the GR relation was applicable to not only ordinary earthquakes 55,58,59 but also LFEs. In this section, the word "earthquake" includes LFEs. We employed the Entire-Magnitude-Range (EMR) technique 60 , which simultaneously calculates the a-and b-values and the completeness magnitude M c , above which all events are considered to be detected. We always paid attention to a robust estimate of M c , because it is critical for analyses of the ETAS models (details of catalog quality evaluation in "Methods"). EMR applies the maximum-likelihood method when computing the b-value to events with magnitudes above M c . Uncertainty in b was according to Shi and Bolt 61 . Substitution of M c , N at M c , and the maximum-likelihood b-value into M, N, and b of a = log 10 N + bM, respectively, gives the a-value. Supplementary Fig. S4 shows a good fit of the GR relation to observations in the present cases.
To compute M c , the EMR technique 60 combines the GR relation with a detection rate function. Details are provided next. Statistical modeling was performed separately for completely detected and incompletely detected parts of the frequency-magnitude distribution. The b-and a-values in the GR relation were computed based on earthquakes above a certain magnitude (M cc ). For earthquakes whose magnitudes were smaller than M cc , it was hypothesized that the detection rate depended on their magnitudes in such a way that large earthquakes were almost entirely detected while smaller ones were detected at lower rates. Earthquakes with M ≥ M cc were assumed to be detected with a detection rate of 1. To evaluate the fitness of the model to data, the log-likelihood was computed by changing the value of M cc . The best fitting model was that which maximized the log-likelihood value.
The software package ZMAP 62 was used to facilitate the computation of a, b, and M c based on the EMR method. Although the package, whose code is open, is written in Mathworks' commercial software language Matlab®, no knowledge is needed since ZMAP is GUI-driven. ZMAP combines many standard seismological tools. A user can use ZMAP to create a graph of frequency-magnitude distribution with the GR relation with a, b, and M c values calculated by EMR ( Supplementary Fig. S4).
ETAS model. The ETAS model 30 was originally introduced for ordinary earthquakes, but we assumed that the model can be extended and applied for LFEs beneath Mount Fuji. In this section, when the word "earthquake" is used, the reader should understand that it also includes LFEs.
The ETAS model is a point-process model that represents the activity of earthquakes of a minimum magnitude (M th ) and above in a certain region during a specified time interval. Seismic activity includes the background activity at a constant occurrence rate μ (Poisson process). The model assumes that each earthquake (including the aftershock of another earthquake) is followed by aftershocks. Aftershock activity is represented by the Omori-Utsu formula 63 in the time domain. The rate of an aftershock occurrence at time t following the i-th earthquake (time t i and magnitude M i ) is given by ν i (t) = K 0 exp{α(M i -M th )}(t-t i + c) -p for t > t i , where K 0 , α, c, and p are constants, which are common to each target aftershock sequence in a region. The rate of occurrence of the whole earthquake series at t becomes (t|H t ) = µ + S<t i <t ν i (t) . The summation is performed for all i satisfying t i < t. Here, H t represents the history of occurrence times with associated magnitudes from the data {(t i , M i )} before time t. The parameter set θ = (μ, K 0 , α, c, p) represents the characteristics of seismic activity. The units of the parameters are day −1 , day −1 , no unit, day, and no unit, respectively. For the case of K 0 = 0, the ETAS model reduces to the Poisson model ( Supplementary Fig. S6). We estimated these parameters using the maximum likelihood method. Because K 0 depends on M in this model, it is necessary to assume a magnitude at which a value for K 0 needs to be known. Throughout this study, M = 2 was assumed for estimating K 0 .
Using the maximum likelihood estimate, it is possible to visualize how well or poorly the model fits an earthquake sequence by comparing the cumulative number of earthquakes with the rate calculated by the model. If the model presents a good approximation of observed seismicity, an overlap with each other is expected. Ordinary time can be converted to transformed time in such a way that the transformed sequence follows the Poisson process (uniformly distributed occurrence times) with unit intensity (occurrence rate) so that visualization can be achieved in two ways [31][32][33] : one graph using ordinary time and the other using transformed time (Fig. 3). Included in the latter graph is the parabola of 95% significance. When the number of earthquakes is insufficiently large, significance actually depends on sample size due to the estimation accuracy of the parameters. The significance of deviation is defined in the case where the empirical curve deviates outside the parabola.
A FORTRAN program package (SASeis2006) associated with a manual for the ETAS analysis was used to calculate maximum likelihood estimates and also to visualize model performance 64 . This has been extended to the program package XETAS 31 using GUI.
When the stationary ETAS model does not fit a dataset well, the simplest alternative model is a "two-stage ETAS model" that considers different parameter values in subperiods before and after a particular time, referred to as change-point T c . AIC is used to test whether or not the changes in seismicity pattern at T c reduces model selection 34 . In this procedure, we separately fitted the ETAS models for each divided period and then compared their total goodness-of-fit values against the one-fit value over the whole period using the principle of minimum AIC. AIC was calculated from the maximum log-likelihood and number of adjusted parameters.
If T c is hypothetically prefixed based on some information other than the occurrence data themselves, such as a notable geophysical event or a notable outside large earthquake, AIC single (AIC for the model fitted over the whole period) can be compared with AIC 2stage (AIC for the 2-stage model fitted on divided periods) to select the model with the smaller value that performs a better fit to the data in the entire target period. If T c is searched from the target data, the 2-stage model becomes harder to accept. Namely, AIC 2stage is compared with AIC single plus the penalty term 2q to select the model. Here q is the degree of freedom to search for the best candidate www.nature.com/scientificreports/ T c from the data. q depends on sample size (number of earthquakes in the target period) 33,65 : q increases with sample size and, for example, lies in the 4-5 range for q between 100 and 1000. Although we adopted AIC for model selection, it cannot always be used for other cases such as the identification of possible earthquake precursors in ionospheric electric content (TEC) 66 .
Time-dependent μ and K 0 . The standard stationary ETAS model can be temporally extended to the nonstationary ETAS model 33,[67][68][69] in such a way that μ and K 0 are assigned as a function of t. The function μ(t) and K 0 (t) are represented by a broken line connecting the respective sequences (t i , μ(t i )) and (t i , K 0 (t i )) for the i-th earthquake, using a Bayesian function.
Although such a sophisticated model is available, a simpler approach was taken to capture essential aspects of the time-dependent background and aftershock activities ( Fig. 4 and Supplementary Fig. S7 and Tables S1-S3). This involves taking a time-window approach. In Fig. 4, the time-windows 2003-2010 and 2012-2019 were considered where the time-window of 2011, which included the Tohoku and Shizuoka earthquakes, was not considered. When creating Fig. 4 (Supplementary Table S1), μ and K 0 were computed, given that other parameters (α, c, p) were constants, where the values for (α, c, p) were obtained as follows. The parameters θ = (μ, K 0 , α, c, p) were first computed for each time-window, each M th , and each CC th . Typical values for α, c, p were α = 0.5, c = 0.0015, and p = 2.8. Using these typical values, parameters μ and K 0 were recomputed for each time-window, each M th , and each CC th . The error bars for μ and K 0 , which can be calculated by XETAS 31 , are based on error distribution depending on the sample size when the number of LFEs is not large enough 65 .
K 0 was forced to express time-dependent aftershock activity (Fig. 4), so the same analysis can be performed for other parameters to support the assumption that the μ-CC th pattern generally remains stable. Namely, p was considered as a time-variable parameter, given that other aftershock parameters (K 0 , α, c) were prefixed as constant, resulting in Supplementary Fig. S7. The same analysis was performed by assuming time-variable α, given that K 0 , c, p were constant ( Supplementary Fig. S7). The parameter values are summarized in Supplementary Tables S1-S3.
The slope (g) and intersection (h) of the least-square regression line and the square of the sample correlation coefficient (R 2 ) for

Data availability
The datasets used and/or analyzed during the current study are available from the corresponding author on reasonable request. The JMA catalog was obtained from https:// www. data. jma. go. jp/ eqev/ data/ bulle tin/ hypo. html. The waveform records were obtained from the permanent stations of the National Research Institute for Earth Science and Disaster Resilience, the Earthquake Research Institute at the University of Tokyo, JMA, and the Hot Springs Research Institute of Kanagawa Prefectural Government. Locations of active volcanoes used for Fig. 1a were obtained from https:// www. mri-jma. go. jp/ Dep/ sei/ fhiro se/ plate/ en. Plate Data. html. The fault model of the 2011 Tohoku earthquake, used to create Fig. 1a, was obtained from Asano et al. 70 . The fault model of the 2011 Shizuoka earthquake, used to create Fig. 1b and Supplementary Fig. S1, was obtained from Fujita et al. 19 . The seismicity analysis software ZMAP 62 , used for Supplementary Fig. S4, was obtained from http:// www. seismo. ethz. ch/ en/ resea rch-and-teach ing/ produ cts-softw are/ softw are/ ZMAP. The program XETAS 31 , used for Figs. 2-4 and Supplementary Figs. S5-S7 and Tables S2 and S3, was obtained from http:// evrrss. eri.u-tokyo. ac. jp/ softw are/ xetas/ index. html. Generic Mapping Tools (GMT) 71 , used for Fig. 1a,b and Supplementary Fig. S1, is an open-source collection (https:// www. gener ic-mappi ng-tools. org).